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A neutral structure in the DD* system around the DD* mass threshold is observed with a 
statistical significance greater than lOcr in the processes e'*'e“ — >■ + c.c. and e^e“ — >■ 

_D°_D*°7r° + c.c. at ^/s = 4.226 and 4.257 GeV in the BESIII experiment. The structure is denoted 
as Zc(3885)°. Assuming the presence of a resonance, its pole mass and width are determined to be 
(3885.7lg y(stat)±8.4(syst)) MeV/c^ and (35lJ^2(®t^t)±15(syst)) MeV, respectively. The Born cross 
sections are measured to be a{e'^e~ —>■ 2t;(3885)°7r°, Zc(3885)° —>■ DD*) = (77± 13(stat)±17(syst)) 
pb at 4.226 GeV and (47 ± 9(stat)±10(syst)) pb at 4.257 GeV. The ratio of decay rates 

0D*0+c determined to be 0.96±0.18(stat)±0.12(syst), consistent with no isospin 

violation in the process Zt;(3885)° —>■ DD*. 

PACS numbers: 14.40.Rt, 13.25.Gv, 13.66.Be 


The existence of exotic states beyond those of con¬ 
ventional mesons and baryons was debated for decades, 
mostly because no convincing experimental evidence 
for them had been found [l|. In recent years, the 
discovery of charged Zc charmonium-like states [1, Q, 
which decay to a charmonium state plus a pion or a 
pair of charmed mesons and, therefore, must consist of 
at least a four constituent quark configuration ccqq', has 
stirred excitement about these possible exotic states. In 
e"'"e“ —>■ TT^ processes, four states have been 
discovered in the decays of Zc(3885)^ —>■ (DD*)^ d, 3, 
Zc(3900)± ^ ^020)± ^ 7r±hc ||, 

and Zc(4025)^ —>■ {D*D*)^ [l^. There have been 
many theoretical predictions and interpretations Q to 
explain their nature as exotic mesons. However, none of 
these models have either been ruled out or established 
experimentally. 

After the discoveries of the charged Z^ states, BESIII 
reported studies of their neutral partners in the isospin 
symmetric channel of e+e“ —>■ tt^Z^. A Zc(3900)° 
is found in e+e“ — TT^-K^J/tp [llj, a Zc(4020)° in 
e+e” —?► TT^TT^hc [m, and a Zc(4025)° in e+e” —>■ 
tt°{D*D*)° [i3|. Evidence for Zc(3900)° in e+e” 

was previously reported with CLEO-c data at 
y/s = 4.17 GeV Q. These measurements indicate that 
the Zc(3900), Zc(4020) and Zc(4025) are three different 
isospin triplet states, since their relative Born cross 
sections of the charged modes to the neutral modes are 
compatible with isospin conservation. This motivates 
a search for the neutral partner of the Zc(3885)^ in 
e"'"e“ {DD*)^^^ + c.c. to identify its isospin. 

In this Letter, the process e+e“ —>■ + c.c. 

is studied, where {DD*)° refers to D^D*~ or D°D*°. 
A neutral charmonium-like structure, the Zc(3885)*^, 
is observed around the (DD*)^ mass threshold in the 
(DD*)^ mass spectrum. This analysis is based on data 
samples collected by the BESIII detector with integrated 
luminosities of 1092 pb“^ at ■>/s = 4.226 GeV and 826 
pb-i at Vs =4.257 GeV [H, [H. Note that charge 
conjugation is always implied, unless explicitly stated. 

BESIII [l^ is a general-purpose detector at the double¬ 
ring e+e“ collider BEPCII, which is used for the study 
of physics in the r-charm energy region [i3- . Monte 
Carlo (MG) simulations based on Geant 4 [l^ are 


implemented in the BESIII experiment. Eor each energy 
point, we generate a signal MG sample based on the 
covariant tensor amplitude formalism [l9j to simulate the 
S-wave process e+e“ —>■ ^ {DD*)^Tr^, assuming 

that the Z^P has = I’*'. Effects of initial state 
radiation are taken into account with the MG event 
generator kkmc [2^ 211, where the line shape of the 
Born cross section of e+e“ —>■ Z^ir^ —>■ {DD*)^Tr^ is 
assumed to follow that of the charged channel e“''e“ —>■ 
Z^TT^ {DD*)^tt^ [3]. In addition, a large statistics 
MG sample of the three body process e+e“ —^ {DD*)^tt^ 
is generated according to phase space (PHSP). To study 
possible backgrounds, MG simulations of V(4260) generic 
decays, initial state radiation production of the vector 
charmonium states, charmed meson production, and the 
continuum process e+e“ qq {q = u, d, s) equivalent 
to 10 times the luminosity of the data at y/s = 4.226 and 
4.257 GeV are generated. Particle decays are simulated 
with EVTGEN pllii for the known decay modes with 
branching fractions set to the world average [ij and with 
the LUNDCHARM model 2^ for the remaining unknown 
decays. 


In this work, we study e+e —>■ D'^D* 7r°, 
D*~ D^tt~ based on the detection of the D'^D^ 
pair and e+e" D°Z)*°7r°, D*° l)°7r° based on 

the detection of the pair. The DD meson pairs 

are reconstructed through five hadronic decay modes 
A'“7r“'"7r“*', Ar“7r"''7r“*'7r°, Ars7r“*", iGsTT+TT*^, K for 
the D^ and three modes AT+tt”, Ar+7r+7r+7r“ 

for the D^. The primary tt^, which is produced along 
with the DD* in the e+e” reaction, is reconstructed 
from a pair of photons, while the soft tt from the D* 
decay is not required to improve the detection efficiency. 
The D'^D~ mode is not included because of its low rate 
compared to D^D^ and D^D^. 

In this analysis, the selection criteria in Ref. @ 
are used to identify the photon, 7r° and ATg 

candidates. The charged-particle tracks in each D 
candidate are constrained to a common vertex, except 
for those from ATg decays, and the of the vertex 
fit is required to be less than 100. Each D candidate 
is required to have its reconstructed invariant mass in 
the range (1.840, 1.880) GeV/c^. Furthermore, a mass- 
constrained kinematic fit (KF) to the nominal D mass is 
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FIG. 1. Distributions of RM(D7r°) at y/s — 4.257 GeV. 
The signal and PHSP processes are overlaid with an arbitrary 
scale. The solid arrows indicate the selection criteria for the 
{DD*)°-k° candidates. Data at y's = 4.226 GeV show similar 
distributions and are omitted. 


performed, and the KF chisquare x% is required to be less 
than 100. In case there is more than one DD combination 
in an event, only the candidate with the minimum sum 
of Xd + Xf) is kept. The DD four-momenta from the 
mass-constrained KF are used for the further analysis. 

The primary tt^ candidates are reconstructed with 
pairs of photons which are not used in forming the DD 
mesons, and their invariant masses M{xi) must be in the 
range (0.120, 0.150) GeV/c^. To reduce backgrounds and 
to improve the resolution, a KF with 2 degrees of freedom 
(2C) is performed, constraining M{xi) to the nominal 7r° 
mass m(7r°) and the recoil mass of ir'^DD, RM{7 t°DD), to 
the nominal tt mass. The 2C KF chisquare X 2 c('^) must 
be less than 200. For each DD mode, if there is more than 
one primary candidate, the one with the minimum 
X 2 c(’’’) retained for further analysis. For e+e” —>■ 
with D*^ —^ Zi°7r°, the process e+e“ —>• 
£) 0 ^*o^o jg major background. To 

reject this background, we require X 2 c('^°) ^ 
perform a similar 2C KF but constrain to 

be zero, which corresponds to the mass of the photon 
in D*° —?► and the corresponding ht chisquare is 

required to satisfy X 2 c( 7 ) > ^0 to further suppress this 
background. The fitted four-momentum of the primary 
7r° is used in the next stage of the analysis. 

In the surviving events, the occurrence of multiple 
{DD*)^tt'^ combinations per event is negligible. To help 
separate the signal events, we require M{D^'k'^) > 2.1 
GeV/c^ and M(Zl°7r°) >2.1 GeV/c^ [2^. Because of the 
limited phase space, the invariant mass of Zl+7r°(Zl°7r°) 
and that of Z)°7r° are highly correlated, and the back¬ 
ground with the selected 7r° and D^ from the D*^ decay 
is suppressed by the above selection criteria, too. The 
RM(Il7r'’) distributions are illustrated in Fig. [1] where 
clear peaks are seen over simulated backgrounds around 
the m{D*) position. These peaks correspond to the 
final states of {DD*)^tt'^. We further require events to 
be within the mass window |RM(Il7r°) — 'm{D*)\ < 36 
MeV/c^ for the final analysis. 


The M{DD*) distribution of the surviving events is 
plotted in Fig. [2] An enhancement near the DD* mass 
threshold around 3.9 GeV/c^ is visible, which is seen in 
both D+D*-7t° and D°D*°7t° at = 4.226 and 4.257 
GeV. As verihed in MC simulations, these structures 
cannot be attributed to the e+e“ —>■ {DD*Y'tt'^ three 
body PHSP or inclusive MC background. Possible 
backgrounds from e+e” —>■ D^*^D** —>■ DD*tt have 
been studied. Most of them, such as D*D*(2'^00), 
DD* (24,60) and D*D* (2420) cannot contribute to the se¬ 
lected events since their mass thresholds are higher than 
4.26 GeV/c^. The only possible peaking background 
e+e“ —>■ D^*^Di(2420) has been studied in Ref. jHl, and 
its contribution is found to be negligible. 

Assuming that there is a resonant structure close to the 
DD* mass threshold (labeled as .Zc(3885)°), we model 
its line shape using a relativistic 5'-wave Breit-Wigner 
function with a mass-dependent width multiplied with a 
phase space factor q 


_y/MFHMj/c^_ 

A/2 -m? + iM(Ti(M) + r2(M))/c2 


a =1,2), 


where Ti(M) = F/ • (m/M) ■ (pj/p'j). I denotes the 
different decay modes, where / = 1 represents the 
D+D*~ decay mode and 1 = 2 represents the D^D*^ 
decay mode. M is the reconstructed mass, m is the 
nominal resonance mass and F/ is the partial width of 
the decay channel I. Under the assumption of isospin 
symmetry, we take F/ to be half of the full width F, 
assuming that the decay rates to other possible coupled 
channels are negligible. p*j(q) is the momentum of the 
D('k^) in the rest frame of the DD* system (the initial 
e“'"e“ system), and pj i® fhe momentum of the D in the 
resonance rest frame at M = m. 

An unbinned maximum likelihood fit is performed on 
the M(DD*) spectra for —?► (/?Z)*)°7r° simultane¬ 

ously at y/s = 4.226 and 4.257 GeV. Three components 
are included in the fits: the Zc(3885)° signal, the 
PHSP processes and MC simulated backgrounds. The 
signal shape is described as a mass-dependent-efficiency 
weighted Breit-Wigner function, described above, con¬ 
voluted with the experimental resolution function. The 
resolution function and the efficiency shape are ob¬ 
tained from MC simulations. The shape of the PHSP 
processes is derived from MC simulations, and their 
amplitudes are allowed to vary in the fits. The inclusive 
MC background distributions are modeled based on 
the kernel estimation , and their sizes are fixed 
according to the expected numbers estimated in the 
inclusive MC samples. The simulated backgrounds are 
validated by comparing their M(DTr^) and RM(/l7r°) 
distributions with those for data in sideband regions 
(1.920,1.974)U(2.090, 2.180) GeV/c^ for the D+D° mode 
and (1.920,1.971) U (2.090, 2.160) GeV/c^ for the D°D° 
mode. 
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FIG. 2. (Upper) Projections of the simultaneous fit to the 
M{DD*) spectra for e''"e“ —>■ D'^D*~'k^ and D^D*°n^ at ^/s 
= 4.226 and 4.257 GeV. (Lower) Sum of the simultaneous 
fit to the M{DD*) spectra for different decay modes at the 
different energy points above. 


We define the ratio TZ = Bd+d— where 

Bd+d*-{Bd°d-°) is the branching ratio of Zc(3885)° —>■ 
D+D*~ D*^). In the fit, TZ is assumed to be 
same for the data at = 4.226 and 4.257 GeV. The 
number of observed signal events, fVobs, is given by 
iVobs = £cr£,£,. (I+(5"'‘*'‘^)(l + 5™‘=)eSint, where cr£,£,. is the 
Born cross section cr(e+e“ —>■ Zc(3885)°7r°, Zc(3885)° —>■ 
DD*), C is the integrated luminosity, (I+J*'^'^) is the 
initial radiative correction factor, (l+J™'^) is the vacuum 
polarization factor 27|, e is the detection efficiency and 
Sint is the product of the decay rates of the intermediate 
states. 

Figure [2] shows the fit results. To assess the goodness 
of fit, we bin the data set in 19 bins such that each bin 
contains at least 10 events, and compute the between 
the binned data and the projection of the fit. We find 
X^/d.o.f. = 18.5/19 for the simultaneous fit in the lower 
plot. The statistical significance of the 2'c(3885)° signal is 
estimated to be more than 12(T, based on the difference of 
the maximized likelihoods between the fit with and with¬ 
out including the signal component. The mass and width 
of the Zc(3885)° are measured to be m(Zc(3885)°) = 
(3894.7 ± 3.0) MeV/c2 and r(Zc(3885)°) = (36 ± 17) 
MeV. The corresponding pole mass and width are calcu¬ 
lated to be mpoie(^c(3885)°) = 3885.7t5;7 MeV/c^ and 


TABLE 1. Summary of systematic uncertainties for the 
resonance parameters, the Born cross sections and the ratio 
of decay rates. Values outside the parenthesis represents 
uncertainties for at ^/s = 4.226 GeV, while those inside 

are for (Jod* a^t ^/s = 4.257 GeV. The total systematic 
uncertainties are obtained by combining all the independent 
sources in quadrature. 


Source 

mpole(MeV/±) 

rpoie(MeV) 

(%) 

77(%) 

Beam energy 

1.0 

3.0 

4 (5) 

1 

Signal shape 

3.5 

8.2 

5 (4) 

2 

Background 

6.8 

6.6 

15 (15) 

4 

Fit range 

0.3 

0.3 

3 (1) 

1 

Mass shift 

3.0 




Resolution 


9.5 

11 (4) 

1 

Efficiency 



11 (11) 

11 

Input-output check 

1.6 

2.5 



(1^5rad)(i_p^vac) 



5 (5) 


Bint 



5 (5) 

5 

C 



1 (1) 


Total 

8.4 

15 

23 (21) 

13 


rpoie(^c(3885)°) = 35ti2 MeV Q. From the fit, we 
determine (Jdd- to be (77 ± 13) pb and (47 ± 9) pb at 
^/s = 4.226 and 4.257 GeV, respectively. We also obtain 
7^ = 0.96 ±0.18. 

The systematic uncertainties on the measurements of 
the Zc(3885)° resonance parameters, the cross section 
(Tdd* and the ratio TZ are studied, and the major 
contributions are summarized in Table HI The systematic 
uncertainties on the Zc(3885)° resonance parameters 
mainly come from the signal shape, background, mass 
shift and detector resolution. The dominant systematic 
uncertainties on and TZ are from the background, 

resolution and detection efficiency. 

The uncertainty from the beam energy is estimated 
by varying the beam energy by ±1 MeV in the 2C 
KF, and the maximum differences of the mass, width, 
at =4.226 (4.257) GeV and TZ are found to 
be 1.0 MeV/c^, 3.0 MeV, 5%(4%) and 1%, respectively. 
To assess the uncertainty of the signal shape, an S- 
wave relativistic Breit-Wigner function with constant 
width [ 2 ^ is taken as an alternative signal model in 
the simultaneous fit. The changes of the fitted mass 
and width are determined to be 3.5 MeV/c^ and 8.2 
MeV, while the change on is 5%(4%) at ^/s =4.226 
(4.257) GeV and on TZ 2%. The systematic uncertainty 
due to the background description is estimated by leaving 
free the absolute numbers of the inclusive backgrounds in 
the fit, or adjusting their shapes by varying the scalings 
of different background components in the inclusive MG 
samples. Those fit results differ from the nominal results 
by 6.8 MeV/c^ in mass, 6.6 MeV in width, 15% in 
(Jqq, both at =4.226 and 4.257 GeV, and 4% in TZ. 
Maximum fluctuations due to changing the ht range are 
assigned as systematic uncertainties. The MG simulation 
of the mass shift and resolution may not fully reflect 
the effects in data, and it is studied by fitting the 
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D* peak in the RM(£)7r°) spectra to obtain the mass 
shift and the resolution difference between data and MC 
simulations. The obtained mass shift is quoted as part of 
the systematic uncertainties of the mass. The variations 
of the fit results after considering the resolution difference 
is assigned as systematic uncertainty. 

Efficiency-related systematic uncertainties are univer¬ 
sal in each D decay mode and include six sources: track¬ 
ing efficiency, particle identification, photon detection 
efficiency, 7r° reconstruction efficiency, Ks reconstruc¬ 
tion efficiency and KF efficiency. The uncertainties 
of tracking efficiency and particle identification for 
and are evaluated to be 1% per track m, 113 . 
The uncertainty in the photon-reconstruction efficiency 
is estimated to be about 1% per photon [3l|. The 
efficiency difference of reconstructing the Ks in MC 
simulations and in data is 4.0% [s^. The uncertainty 
in 7r° reconstruction is 1% jsij. The systematic bias 
of the KF is estimated by using the track-parameter- 
correction method [s^. The correction factors for helix 
track parameters are determined from the control sample 


e^e 


^ K*{SmfK+'K- 


The total 

efficiency-related systematic uncertainty is taken as the 
square root of the quadratic sum of the individual 
uncertainties. The potential bias from the event selection 
and the analysis procedure is studied with input-output 
checks, which compare the output results with the input 
values of the resonance mass and width based on MC 
simulations. We assign the systematic uncertainty of 
1.6 MeV/c^ in mass and 2.5 MeV in width accordingly. 
The systematic uncertainty of the radiative correction 
factor 1 -f (5''*^'^, which includes the effect on the detection 
efficiency, is estimated to be 5% by changing the input 
(ZlZ)*)°7r° line shape within errors [^. The systematic 
uncertainty of the vacuum polarization factor 1 -|- is 
0.5% taken from the QED calculation 27|. The weighted 
systematic uncertainty of Sint is from the world average 
value [l|. The uncertainty of integrated luminosity is 
taken as 1% by measuring Bhabha events [ 3 . The 
uncertainty of the mass window requirement is negligible. 
The overall systematic uncertainties are determined by 
combining all the sources in quadrature, assuming they 
are independent. 

In summary, we study e“'"e“ —>■ D'^D*~'it^ + c.c. 
and e+e" —>■ + c.c. using data taken at 

y/s = 4.226 and 4.257 GeV. A neutral structure 
around the DD* mass threshold is observed with a 
statistical significance greater than lOcr. Assuming 
that it is a resonance, we model it with a relativistic 
Breit-Wigner function. Its pole mass and width are 
measured to be (3885.7l5'7(stat)±8.4(syst)) MeV/c^ 
and (35;i2 (stat)±15(syst)) MeV, respectively, which 
are close to the mass and width of the reported 
charged Zc(3885)+ 0 i- The Born cross sections 


cr(e 


+ c“ 


—>■ {DD*y-K^ + c.c.) are determined to 


be (77 ± 13 ± 17) pb and (47 ± 9 ± 10) pb at = 4.226 


and 4.257 GeV, respectively, which are consistent with 
half of (7(e+e“ —>■ — ?► {DD*)^'k~ -|- c.c.) 0|. A 

comparison between the resonance parameters of the 
Zc(3885)“'' and the ^^(3885)° is summarized in the 
Supplemental Material [ 2 ^ . All these observations 
favor the assumption that the ^c(3885)° is the 
neutral isospin partner of the Zc(3885)^, and the 
Zc(3885)^/^c(3885)° form an isospin triplet. In 
addition, we determine the ratio of the decay rate 

= 0.96 ± 0.18 ± 0.12, which is 


^ _ B(Ze(3885)°->£)+£>*-) 


B(Zc(3885)«->D“D*0) 

consistent with unity. Hence, no isospin violation in the 
process Zc(3885)° —^ DD* is observed. 
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